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CHAPTER  1 
INTRODUCTION 


In  recent  years  much  of  the  work  done  in  ccmputational  fluid  dynamics  (CFD) 
has  been  to  produce  high  resolution  flow  solvers  !25.  31,  32  and  to  develop  domain 
decomposition  techniques  m  order  to  accurately  model  the  aerodynamics  of  com.plex 
geometries  4,  5,  11,  12,  13,  24;,  Although  these  high  resolution  sol-ers  can  pro¬ 
vide  good  answers  they  are  still  grid  dej  indent  and  will  usually  require  significant 
computer  time  to  converge  for  all  but  the  simplest  problems. 

Grid  dependency  can  be  taken  care  of  by  generating  “good”  grids,  where  “good  ’ 
grid  simply  means  any  grid  that  causes  the  numerical  solution  to  converge  rapidly 
to  the  actual  physical  solution  for  a  given  geometry.  Since  most  real  geometries  are 
not  necessarily  simple  it  is  difficult,  if  not  impossible,  to  generate  a  single  “good” 
grid.  Several  techniques  have  been  developed  to  allow  the  domain  of  interest  to  be 
approximated  by  a  combination  of  smaller  “good”  grids  as  opposed  to  using  a  single 
complex  grid.  These  domain  decomposition  techniques  have  been  applied  to  the 
analytic  solution  of  partial  differential  equations  for  many  years. 

One  of  the  simpler  forms  of  domain  decomposition  is  referred  to  as  composite 
blocked  or  multiblocked  grid  structures  |4,  28,  29,  30|.  In  this  method  each  block 
(separate  grid)  covers  a  given  region  of  the  domain  and  the  interfaces  between  these 
blocked  regions  match  physical  grid  points  exactly,  although  the  slopes  of  the  grid  lines 
through  these  boundaries  may  be  discontinuous.  In  fact  when  a  singularity  is  present 
on  a  block  boundary  multiple  lines  from  one  block  may  merge  at  the  singularity  and 
become  a  single  grid  line  in  the  adjoining  block. 

Another  widely  used  method  is  patched  or  zonal  grids  '18,  23,  24].  For  zona!  grids 
the  boundaries  between  the  zones  must  coincide;  however,  the  grid  lines  i  tersecting 
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these  boundaries  do  not  cross  the  boundary  in  a  continuous  manner.  In  fact  the 
common  boundary  between  zones  may  not  contain  the  same  number  of  points.  One 
imnortant  restriction  with  this  grid  system  is  that  a  sufficient  number  of  points  must 
be  on  the  zonal  boundaries  to  maintain  conservation  principles  across  such  boundaries. 
Blocked  grids  are  a  special  case  of  zonal  grids. 

.■\n  even  more  general  domain  decomposition  method  is  referred  to  as  embedded, 
overset,  or  Chimera  grids  13,  5.  13.  14,  22:.  When  using  embedded  grids  no  set  bound¬ 
ary  correspondence  ia  required  betw'een  grids.  Grid  boundaries  that  lie  in  the  interior 
of  the  domain  are  set  by  interpolation.  This  method  also  requi’^es  sufficient  grid  points 
for  intergrid  communication  to  occur  in  a  conservative  manner.  Grid  points  that  are 
within  the  overlapped  region  and  points  that  fall  inside  of  solid  surfaces  are  not  used 
as  a  part  of  the  solution.  These  points  are  referred  to  as  blanked  or  1- blanked.  Since 
these  grids  are  generated  independently  they  can  be  added,  deleted,  or  moved  in 
order  to  make  configuration  changes  without  having  to  regrid  the  entire  domain,  A 
great  deal  of  work  has  gone,  into  automating  the  search  for  the  intergrid  interpolation 
stencils  and  the  grid  movement  procedure  for  grid':  in  relative  motion  13,  22,.  Both 
zonal  and  blocked  grids  are  special  cases  of  embedded  grids. 

These  domain  decomposition  techniques  allow  the  solution  of  more  complex  ge¬ 
ometries:  however,  the  amount  of  computer  timie  required  usually  increases  with  the 
number  of  artificial  boundaries  introduced  into  the  domain.  This  has  lead  to  the 
development  of  several  techniques  for  in.  '  ^ving  convergence  of  the  flow  solvers. 

One  method  for  improving  convergence  of  iterative  schemes  for  solving  systems 
of  equations  is  multigrid  .2,  9,  10,  15,  16.  211.  Regular  multigrid  simply  uses  several 
levels  of  increasingly  coarser  grids  covering  the  same  domain  to  solve  a  set  of  algebraic 
equations.  Tne  solution  is  generated  on  the  finest  level  and  corrections  to  the  solution 
are  calculated  on  the  coarser  levels.  Multigrid  techniques  applied  in  CFD  have  shown 
significant  im.provements  in  convergence  over  single  level  methods  '2,  3,  20.  27.  33j. 
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Multigrid  techniques  have  also  been  developed  to  provide  increased  resolution  in 
portions  of  the  domain.  These  local  refinement  techniques  [9,  11.  15,  21'  usually 
generate  the  local  grid  levels  by  quartering  (in  two  dimensions)  existing  cells  in  the 
global  grid  and  storing  this  refinement  as  a  part  of  the  next  finer  level.  Boundary 
information  on  the  local  level  is  usually  obtained  by  interpolation  from  a  lower  level. 
Accurate  communication  between  these  levels  is  critical  !21j.  These  methods  are 
typically  automated  such  that  additional  levels  are  added  during  the  solution  process. 

Multigrid  has  also  been  used  to  provide  increased  convergence  on  blocked  and  em¬ 
bedded  grid  systems.  When  this  technique  is  applied  to  an  embedded  grid  system, 
typically  a  complete  multigrid  cycle  is  made  in  each  grid  independent  of  the  other 
grids  13)  and  the  communications  between  grids  is  achieved  on  the  finest  grid  level. 
Special  treatment  of  blanked  points  is  required  on  the  the  coarser  levels.  Numerical 
results  show  that  when  multigrid  is  applied  to  the  blocked  grid  system  it  is  better 
to  make  calculations  and  communicate  between  blocks  on  each  level  than  to  perform 
a  full  m.ultigrid  cycle  m  each  block  independently  [331.  This  would  seem  to  indicate 
that  if  multigrid  could  be  applied  in  a  way  that  would  allow  the  low’er  levels  of  the 
embedded  grids  to  communicate  that  even  better  convergence  rates  would  result.  Ref¬ 
erence  17;  applies  multigrid  to  an  embedded  grid  system  to  solve  an  elliptic  problem 
in  this  fashion. 

The  purpose  of  this  study  is  to  show-  the  feasibility  of  using  a  nonaligned  multi¬ 
grid  technique  to  both  communicate  between  embedded  grids  and  to  simultaneously 
increase  the  convergence  rate.  Nonaligned  multigrid  [I5i  can  be  thought  of  as  the 
global  application  of  multigrid  to  a  general  embedded  grid  system.  In  this  case  the 
only  blanked  points  will  lie  inside  of  solid  surfaces  and  all  communication  between 
grids  and  grid  levels  will  be  done  through  the  m.ultigrid  techniques.  This  means  that 
overlapped  regions  will  contain  the  solution  on  one  grid  and  the  other  grids  occupy¬ 
ing  the  same  region  will  be  used  to  generate  the  corrections  to  that  solution.  The 


4 


nonaligned  multigrid  is  basically  a  local  refinement  technique  with  the  local  levels 
generated  independently  of  the  global  level.  This  method  will  be  applied  to  the  ex¬ 
plicit  flux  difference  splitting  [26,  32[  approximation  of  the  two  dimensional  Euler 
equations  for  inviscid  compressible  flow. 

In  Chapter  2  the  nondimensional  Euler  equations  and  the  flux  difference  splitting 
(EDS)  algorithm  [26,  32]  are  presented  in  terms  of  a  body  conforming  curvilinear  co¬ 
ordinate  system.  These  partial  differential  equations  are  stated  in  strong  conservation 
law  form  and  in  quasi-linear  form.  The  latter  of  these  forms  is  necessary  in  order  to 
obtain  the  eigenvalues. and  eigenvectors  of  the  flux  Jacobian  matrices  [4,  31!  that  are 
needed  for  the  EDS  algorithm.  Stability  requirements  for  the  explicit  EDS  ^”d  initial 
conditions  are  stated.  A  phantom  point  formulation  is  also  examined  for  tue  solid 
surface  and  supersonic  inflow /outflow  boundary  conditions. 

Multigrid  methods  are  discussed  in  Chapter  3  and  the  formulation  for  the  Eull 
Approximation  Scheme  (EAS,  sometimes  referred  to  as  the  full  approximation  stor¬ 
age)  [2,  9,  21]  is  presented.  The  use  of  bi-linear  interpolation,  applied  in  the  compu¬ 
tational  domain,  for  the  restriction  and  prolongation  operations  is  examined.  Grid 
scheduling  and  boundary  conditions  for  different  levels  will  also  be  discussed. 

Chapter  4  presents  a  discussion  of  the  embedded  grid  approach.  This  includes  an 
examination  of  the  stencil  jumping,  boundary  interpolation,  and  blanking  procedures. 

In  Chapter  5  the  nonaligned  multigrid  (NAM)  approach  is  described  using  EAS 
for  the  nonlinear  problem.  The  use  of  bi-linear  interpolation,  applied  in  the  physi¬ 
cal  domain,  for  the  restriction  and  prolongation  operations  is  examined.  Grid  level 
ordering  and  boundary  conditions  will  also  be  discussed  for  the  nonaligned  grids. 

A  discussion  of  the  results  for  the  5°  ramp  and  the  5°  ramp  near  a  flat  plate  is 
contained  in  Chapter  6.  These  include  comparisons  between  single  grid,  multigrid, 
embedded  grid,  and  nonaligned  multigrid  calculations.  Comparisons  will  be  made  by 
looking  at  solutions,  convergence  histories,  and  CPU  times  for  various  methods. 
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Conclusions  and  observations  about  the  numerical  techniques  used  in  this  study 
are  presented  in  Chapter  7. 


CHAPTER  2 

EULER  EQUATIONS  AND  FLUX  DIFFERENCE  SPLITTING 


This  chapter  states  the  strong  conservation  law  form  of  the  Euler  equations  in 
two  dimensions  in  terms  of  a  body  conforming  curvilinear  coordinate  system.  The 
equation  is  converted  to  the  quasi-linear  form  and  the  eigenvalues  and  eigenvectors 
are  given  for  the  flux  Jacobian  matrices.  The  first  order  FDS  algorithm  used  to 
approximate  the  Euler  equations  and  the  resulting  stability  condition  are  presented. 
Time  stepping,  initial  conditions,  and  boundary  conditions  used  in  this  study  are  also 
discussed. 


Euler  Equations 


The  equations  of  motion  in  strong  conservation  law  form  for  inviscid  compressible 
flow  are  presented  here  for  two  dimensions.  These  equations  are  nondimensional  and 
are  expressed  in  terms  of  a  general  curvilinear  (body  fitted)  coordinate  system  [4,  26, 
31]. 


dt  dT) 


(1) 


Where  the  primary  variables  Q  and  the  flux  vectors  Fk  [k  is  (  or  rj  denoting  the 
curvilinear  component)  are  given  by 


p 

pOk 

Q  =  J 

pu 

pv 

;Fk  =  J 

pu9k  +  kxp 
pv6k  +  kyp 

e 

.  {e  +  p)6k  . 

In  this  notation  p  is  the  fluid  density,  u  and  v  are  the  Cartesian  components  of 
velocity,  p  is  the  pressure,  and  e  is  the  total  energy  per  unit  volume.  Since  there  are 
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only  four  equations  and  five  unknowns  the  perfect  gas  assumption  is  made  in  order 
to  relate  the  energy  and  pressure. 

e  =  (3) 

The  remainder  of  the  terms  in  Eq.  2  are  the  metrics,  kx  and  ky,  of  the  transformation 
from  Cartesian  to  curvilinear  coordinates,  the  Jacobian  of  this  transformation,  J ,  and 


the  contravariant  velocity  components  6^. 

i.  =  (4) 

=  ~  J  Xrj  (5) 

Vx  =  (6) 

Tjy  =  (7) 

J  =  X(_yr,  -  J/^Ir,  (8) 

6k  =  y-kx  +  vky  (9) 


The  nondimensionalization  [4,  26]  of  Eq.  1  is  obtained  from’ the  following  defi¬ 
nitions.  First,  all  dimensional  quantities  are  indicated  by  a“  ’  ”  and  a  convenient 
length  I  is  defined.  The  quantities  subscripted  by  oo  denote  reference  values  in  the 
undisturbed  gas. 


(10) 


'Where  c^o  =  {iPool PooY^'^ ,  is  the  speed  of  sound  in  the  undisturbed  gas.  Applying 
this  nondimensionalization  to  the  Euler  equations  (Eq.  1)  yields  exactly  the  same 
equation  in  terms  of  nondimensional  quantities. 

The  quasi-linear  form  of  the  Euler  equations  is  obtained  from  Eq.  1  by  defining 
the  flux  Jacobian  matrices  as  Ak  =  dFk/dQ  and  applying  the  chain  rule. 


dQ 

dt 


=  0 


(11) 


When  evaluated  ,4^  becomes 


’  0  fci  ky  0 

kj-cj)  —  uOk  Ok  +  kx(2  —  'y)u  kyU  —  k^i'y  —  l)u  kj.{y  —  1) 

k^jd)  —  vOk  kxV  —  kJj  —  l)u  Ok  T-  ky(2  —  'y)v  kJ^  —  1) 

A4>-^l3)0k  -k,0-u{')-\)Ok  -kytS-vi'-i-DOk  yOk 


with  3  =  (f)  —  7e/p,  and 


3  = 


(7  -  1) 


+  v^] 


(12) 


(13) 


The  eigenvalues  and  the  left  and  right  eigenvectors  of  this  matrix  are  needed  in  order 
to  implement  the  FDS  algorithm  presented  later  in  this  chapter.  In  the  following 
equations  Xk  are  the  eigenvalues,  the  columns  of  Tk  are  the  right  eigenvectors,  and 
the  rows  of  7*.“^  are  the  left  eigenvectors  [4,  31], 


=  Ok 

(14a) 

=  Ok 

(14b) 

=  Ok  +  c{kj  ^  ky^Y^'^ 

(14c) 

K 

=  Ok  -  c{kj  +  ky^Y^^ 

(14d) 

Tk-^ 


1 

■  1  0 

P 

V2c 

p 

V2c 

Tt  - 

U  —pky 

k  — 

V  pkx 

“  ^v) 

.  P{^'^x  -  uky) 

P  (  0+c’ 

C^-<^ 

w(-r-l) 

ri  . 

c2 

uky  —vkx 

_ 

kx. 

T{(f>-0kc]  -  (7  -  l)u)  T(^„c  -  (7  -  l)u) 

T{<f)  +  0kc)  -T(iiC  +  (7  -  l)ii)  -T(^yC  +  (7  -  1)t;) 


k) 


0 

T(7-1) 

T(7-1) 


(15) 


(16) 


Where  indicates  division  by  (k^^  +  ky^Y^^  and 
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First  Order  Flux  Difference  Splitting 


Now  that  the  eigenvectors  and  eigenvalues  of  the  flux  Jacobian  matrix  have  been 
defined  the  flux  difference  splitting  can  be  used  to  approximate  Eq.  1.  The  following 
equation  is  the  explicit  finite  volume  discretization  [25,  26,  32]  used  for  this  study. 


0  (19) 


Where  /j.  is  the  numerical  flux  in  the  k  direction  and  =  Arj  =  1  for  the  uniform 
curvilinear  coordinate  system.  For  a  first  order  approximation  the  numerical  flux  is 


given  bv 


(fiUin.,  =  (/<)m  + 

m—l 

(7,)., ,.1/2  = 


(20a) 


(20b) 


The  flux  vector  F*  is  evaluated  at  (z,  j)  to  obtain  fk  and  A  is  the  non-positive 


eigenvalue  given  by 


The  subscript  m  on  Tk  indicates  the  corresponding  right  eigenvector  contained  in  the 
column  of  the  matrix  defined  in  Eq.  15.  The  a*  term  in  Eq.  20  is  defined  for 
the  first  order  scheme  as 


(22a) 

{a,U  =  -  qZ)i-  (22b) 

In  this  case  the  subscript  m  indicates  the  left  eigenvector  corresponding  to  the  eigen¬ 
value  (and  right  eigenvector)  previously  discussed.  The  left  eigenvector  is  contained 
in  the  row  of  the  matrix  defined  in  Eq.  16.  The  subscript  /  indicates  mul¬ 
tiplication  of  the  row  vector  of  left  eigenvectors  by  the  column  vector  of  differenced 
dependent  variables. 
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These  eigenvalues  and  eigenvectors  given  by  Eqs.  14,  15,  and  16  in  the  previous 
section  are  evaluated  using  the  “Roe”  averaged  variables  defined  by 


P  =  \/PxPx  +  \ 

v/ Pt+l^i+1 
U  =  - — - ; - - 

^JPt  VPi+i 

\/A«i  +  V'Pt  +  ll'.+I 
V  =  - — - ; - - 

x/Pt  +  y/Pt+1 

=  (7  - 

H  =  -(e  +  p) 

p 

p  +  v^) 

e  =  --^+P— T 


(23a) 

(23b) 

(23c) 

(23d) 

(23e) 

(23f) 


These  equations  are  for  the  (  coordinate  direction.  The  terms  are  evaluated  at  the  i 
index  indicated  and  j.  The  equations  for  the  rj  direction  are  obtained  by  substituting 
j  for  i. 

The  “Roe”  averaging  is  selected  to  ensure  uniform  validity  across  discontinu¬ 
ities  [25].  By  using  Roe’s  numerical  fluxes  shocks  and  contact  discontinuities  can 
be  captured  in  as  few  as  one  to  two  grid  cells  [26,  32].  This  method  also  allows  the 
use  of  fewer  grid  points  than  the  flux  vector  splitting  schemes  to  capture  viscous  shear 
layers.  The  two  (and  three)  dimensional  applications  of  this  method  assumes  that 
the  characteristic  waves  propagate  in  a  direction  normal  to  the  grid  interfaces  [32]. 


Choice  of  Time  Step 


Applying  the  von  Neumann  stability  analysis  [1]  to  Eq.  19  yields  the  following 
condition  for  stability  (see  the  Appendix),  analogous  to  the  one  dimensional  Courant 
Friedrichs- Lewy  (CEL)  condition  with  the  Courant  number,  {CFL),  given  by 

CFL  =  A((|«(|  +  +  {J  +  +  >>2)  <  1.  (24) 

Two  variations  on  this  equation  are  implemented.  The  first  is  known  as  global 
time  stepping  and  is  necessary  for  time  accurate  calculations.  In  this  approach  the 
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time  step,  At,  is  chosen  to  satisfy  the  stability  condition  and  all  cells  are  advanced 
the  same  amount  of  time  each  iteration.  In  the  local  time  stepping  formulation  the 
Courant  number  (  labeled  C F L)  is  selected  and  held  constant.  This  method  is  only 
valid  for  steady  state  problems  because  each  cell  may  be  advanced  by  a  different  time 
step. 


Initial  and  Boundary  Conditions 


The  entire  domain  is  initial!}'’  assumed  to  be  at  free  stream  conditions.  The  Mach 
number  is  specified  and  the  remainder  of  the  free  stream  values  are  determined  as 
follows  for  the  nondimensionalized  variables. 


=  ] 

(25a) 

{p'^)oo 

=  Moo 

(25b) 

{pv)oo 

=  0 

(25c) 

^oo 

1 

'>'(7  - 

a. 

-1) ' 

(25d) 

Poo 

=  (7- 

l)(eoc 

-  \mI) 

(25e) 

The  boundary  conditions  are  applied  using  an  explicit  image  point  or  phantom 
point  formulation  [31].  This  formulation  works  by  selecting  a  set  of  fictitious  points 
just  outside  of  the  domain  with  each  cell  center  on  the  boundary  having  a  correspond¬ 
ing  phantom  point.  Cell  centers  at  the  corners  will  have  three  phantom  points,  one 
for  each  coordinate  direction  and  a  corner  point  to  fill  out  the  boundary.  This  corner 
phantom  point  is  not  required  or  used  by  the  FDS  algorithm.  Eqs.  26,  27,  and  28  are 
used  to  set  the  values  of  these  phantom  points  for  the  specified  boundary  conditions. 
In  these  equations  the  subscript  a  denotes  the  phantom  point  and  subscript  b  the 
corresponding  point  in  the  domain. 


12 


Supersonic  inflow 


Pa 

=  1 

(26a) 

{pu)a 

=  A/oo 

(26b) 

iP^)a 

=  0 

(26c) 

Ca 

II 

1 

1 - * 

+ 
to  1 
> 

(26d) 

Pa 

(26e) 

Supersonic  outflow 


Pa 

=  Pb 

(27a) 

{P^)a 

=  {P^)b 

(27b) 

{P^)a 

=  {P^)b 

(27c) 

Sa 

=  Sb 

(27d) 

Pa 

=  Pb- 

(27e) 

Solid  surface  (zero  pressure  gradient) 


Pa  Pb 

(28a) 

{pu)a  =  {pu)b  -2{pek'kx)b 

(28b) 

{pv)a  =  {pv)h  -  2{pekky)b 

(28c) 

Ca 

(28d) 

Pa  =  Pb- 

(28e) 

CHAPTER  3 

MULTIGRID  APPROACH 


The  idea  behind  multigrid  is  to  accelerate  the  convergence  of  a  relaxation  scheme 
by  adding  corrections  to  the  fine  grid  solution  based  on  solutions  generated  on  coarser 
grids  [9,  10,  21j.  In  order  to  describe  why  multigrid  works  two  important  pieces  of 
information  need  to  be  recognized.  First,  most  relaxation  schemes  eliminate  the 
high  frequency  error  in  as  little  as  2  or  3  iterations,  while  the  same  scheme  can  take 
thousands  of  iterations  and  may  never  satisfactorily  eliminate  the  lower  frequenc}’ 
errors.  Secondly  the  low  frequency  errors  on  a  fine  grid  appear  as  high  frequence- 
errors  on  coarser  grids.  It  would  then  appear  that  nearly  all  of  the  error  could  be 
treated  as  high  frequency  error  when  the  correct  sequencing  of  grids  is  used.  The 
following  is  a  description  of  the  multigrid  method  used  in  this  work. 

Full  Approximation  Scheme  (FAS) 


The  full  approximation  scheme  [2,  8,  9,  16]  is  used  a  great  deal  in  computational 
fluid  dynamics  (CFD)  due  to  its  ability  to  handle  nonlinear  problems  such  as  the  Euler 
or  Navier-Stokes  equations.  FAS  restricts  both  the  residual  and  dependent  variables 
to  the  coarser  grids.  This  differs  from  the  multigrid  designed  for  linear  systems  of 
equations  which  only  restricts  the  residual  to  the  coarser  grid. 

Letting  a  symbol  superscripted  by  h  denote  that  the  symbol  belongs  to  the  fine 
grid  level  then  the  grid  sequence  is  denoted  by  superscripts  (/i,  2h,  Ah, .  .  .  ,  nh).  Where 
2h  indicates  a  grid  over  the  same  domain  as  h  with  half  the  number  of  grid  points 
in  each  direction  and  recursively  for  each  subsequent  grid.  In  order  to  apply  FAS, 
Eq.  19  is  written  in  the  following  form 

AQ'*)  =  0  (29) 
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where  is  the  exact  steady  state  solu*^ion  to  the  finite  difference  equation  and  is 
the  finite  difference  operator  on  the  fine  grid,  defined  for  the  explicit  flux  difference 
splitting  algorithm  as 

=  (/f)?.,,,.,  -  (7f)^.,2,,  +  (7,)::, *1,2  -  (7,t-.,2.  (30) 

Applying  an  iterative  scheme 

L\q^)^-R^  (31) 

where  is  the  residual  obtained  by  subtracting  Eq.  31  from  Eq.  29. 


L^iQ^)  -  L\q^)  =  -R^ 

This  equation  is  then  approximated  on  the  next  coarser  grid  as 


=  IJI\-R^)  -  L^^il^q  )  (33) 

with  the  restriction  operators  11^^  and  to  be  discussed  in  the  next  section  Eq.  31 
is  then  used  to  substitute  for  R^. 

-  Ill\L^q^)  (34) 

The  right  hand  side  of  this  equation  is  now  the  defect  correction  (or  the  relative 

truncation  error)  between  the  solution  on  finest  (level  h)  and  the  coarse  (level  2h) 
grids.  This  defect  correction  acts  as  a  forcing  function  on  the  coarser  levels  in  order  to 
maintain  the  accuracy  of  the  finest  level.  Note  that  if  the  L^[q^)  =  0,  then  q^  = 
and  the  exact  solution  for  in  Eq.  34  is  Designating  the  defect 

correction  as  the  previous  equation  becomes 

=  T^'*.  (35) 

Equations  for  any  number  of  levels  can  be  generated  by  this  same  procedure.  For 
example,  treating  Eq.  35  with  the  same  procedure  used  on  Eq.  29  the  equation  on 
the  Ah  level  becomes 

/.‘•'•(Q*''*)  =  t'*'*.  (36) 
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The  defect  correction  on  this  level  is  actually  a  measure  of  the  relative  truncation 
error  between  level  h  and  ih. 

^  _  112^^  (37) 

Observing  Eqs.  29,  35,  and  36  it  can  be  seen  that  the  same  system  of  equations 
are  solved  on  each  level  (r^  =  0)  making  FAS  an  easy  technique  to  program  (see 
references  ;2,  8]  for  more  details).  These  equations  also  indicate  that  the  residual 
(obtained  from  the  iterative  procedure)  minus  the  defect  correction  is  driven  to  zero 
on  the  coarse  grids  in  the  FAS  procedure.  The  only  real  drawback  to  this  procedure 
is  the  need  to  store  the  solution  at  all  levels.  Fortunately  each  grid  level  requires  only 
d  quarter  (in  two  dimensions)  of  the  previous  level’s  storage  requirement. 

The  remainder  of  the  FAS  procedure  is  to  calculate  the  coarse  grid  correction.  1'. 
and  pass  it  to  the  next  finer  level.  The  following  equations  are  used  to  cilcr:ate  the 
correction  on  the  ‘ih  level  and  to  prolongate  (interpolate)  it  to  the  2h  level. 

(38) 

^  ^2h_^j2hy4h  ^39, 

Restriction  and  Prolongation  Operators 


These  operators  simply  define  the  interpolation  needed  to  move  quantities  between 
grids.  The  type  of  interpolation  to  be  used  is  fairly  arbitrary  and  th^e  is  no  require¬ 
ment  that  the  restriction  and  prolongation  operators  be  of  the  same  type  o:  order.  It 
is  actually  better  to  chose  a  prolongation  operator  that  will  ne  t  pass  high  frequency 
errors  due  to  the  increased  iteration  -  ost  on  the  finer  grid  levels  [9]. 

The  restriction  operator  for  a  point  quantity  given  by, 


1.2J-1 


^2t.2]  - 


+  92t- 


1.2j 


92i,2j  J 


(40) 


is  bi-linear  interpolaaon  applied  in  the  computational  domain.  Since  the  computa¬ 
tional  domain  is  uniform  all  of  the  coefficients  are  equal  and  constant.  Defining  the 
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restriction  operator  in  this  way  will  allow  conservative  transfer  of  mass,  momentum, 
and  energy  only  when  the  grid  is  uniform.  The  subscripts  {i.j)  denote  a  genera!  point 
in  the  coarser  grid. 

The  restriction  operator  for  the  residual, 


R  )ij  —  (•^2t-1.2j-l  ^  ^2i-  l,2j  ' 


(4i: 


sim.pK’  sums  the  residuals  on  the  fine  grid  that  lie  inside  of  thj  coarse  grid  cell.  This 
formulation  gives  an  exact  conservative  transfer  of  residuals  !2:  (reference  2  also  uses 
a  conservative  restriction  operator  for  the  conserved  variables). 

Bi-linear  interpolation  applied  in  the  computational  domain  is  also  used  for  the 
prolongation  operator.  In  this  case  point  quantities  are  transferred  back  to  the  finer 
grid  levels  using 


i^2h^  ^^)2t-\.2]-l 

(4V'^'‘)2.,22 


^(9V;y-3(V 
-  3(V 
]^(9K?  +  3(V 
]4(9V;“  +  3(1 


■2h 

'-1.2 


■2h 

'-1.2 


■2h 

>-1.2 

■2h 

'-1.2 


'i-lj-W 

^  y2/»  \  ^  T  •2>i  S| 

’  '  >,2-1  )  '  '-1,2-1  ) 

^  ^  y2h  \ 

‘  )  *  1-1,2-!  - 

-1.  1  '2*  N  ^  \  ’2h  \ 


(42a) 

,42b) 

(42c) 


(42d) 


Boundary  Conditions 


The  boundary  conditions  applied  to  the  fine  grid  are  given  in  Chapter  2,  with  the 
exception  of  the  corner  phantom  points  which  are  set  by  taking  the  average  of  the 
two  nearest  phantom  points.  These  points  are  not  used  in  the  flux  difference  splitting 
but  they  do  come  into  play  during  the  transfer  of  information  between  grid  levels. 

The  fine  grid  boundary  conditions  cannot  be  applied  directly  to  the  coarser  grid 
levels  without  causing  large  truncation  errors  [16].  This  can  be  corrected  by  a  bound¬ 
ary  condition  defect  correction  generated  in  a  manner  similar  to  the  defect  correction 
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for  the  interior  region.  The  problems  with  this  are  that  the  boundary  conditions  must 
have  a  flux  formulation,  the  calculation  of  the  correction  itself  can  be  very  complex, 
and  a  correction  of  this  type  can  make  the  program  problem  dependent  16,.  An 
alternative  to  this  is  to  pass  the  boundary  conditions  from  the  fine  grid  then  perform 
relaxation  on  the  boundaries  of  the  coarser  levels  !9|.  Another  option  is  to  update  the 
boundary  conditions  only  on  the  fine  grid  and  freeze  them  on  the  coarser  levels  [161. 
For  finite  volume  schemes  using  a  phantom  point  boundary  condition  formulation 
simpl}’  applying  the  fine  grid  conditions  to  the  coarse  grid  is  acceptable  when  only 
the  change  on  the  coarse  grid  is  transferred  to  the  fine  grid  [16].  This  last  method 
of  applying  boundary  conditions  was  .selected  in  order  to  avoid  the  need  for  special 
interpolation  (and  extrapolation)  that  would  be  necessary  to  update  the  phantom 
points. 


V- Cycle  Procedure 


Level 

S)  (T)  h 

2h 

4h 

8h 

Figure  1:  FAS  V-Cycle  Grid  Schedule 

Several  methods  of  scheduling  levels  can  be  used  for  multigrid  [2,  10].  The  \'-Cycle 
procedure  is  used  in  this  study.  This  procedure  is  recursive  and  is  easily  applied  to 
any  number  of  levels.  The  actual  grid  scheduling  for  a  four  level  scheme  is  shown  in 
Figure  1  and  the  corresponding  solution  procedure  follows. 
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1.  Iterate  Vh  times  on  =  —At{L^(q^)) 
for  (  n  =  2h;  n  <  8h;  n  =  2*n  )  { 

2.  Calculate  the  residual 

3.  Restrict  solution  to  coarser  level  g”  =  ■^n/29”^^ 

4.  Restrict  residual  minus  defect  correction  to  coarser  level  ~ 

(recall  =  0) 

5.  Calculate  the  defect  correction  t"  =  -^"(-^n/29”^^)  ~ 

6.  Iterate  times  on  =  —At{L^{q^)  —  r”) 

} 

for  (  n  =  8h;  n  >  2h;  n  =  n/2  )  { 

7.  Calculate  the  coarse  grid  correction(CGC)  ]/"  =  Q"  — 

8.  Prolongate  CGC  and  add  to  solution  on  finer  level 

9.  Iterate  i/„  times  on  Ag"^^  =  —  —  r"^^) 

} 


10.  If  solution  is  not  converged  go  to  1 


CHAPTER  4 

EMBEDDED  GRID  APPROACH 


Embedded  grids  are  also  referred  to  as  Chimera  or  overset  grids.  This  domain  de¬ 
composition  method  uses  independently  generated  grids  to  accurately  model  different 
regions  of  the  domain.  Since  these  grids  are  generated  independently,  configuration 
modifications  will  not  require  regridding  the  entire  problem.  This  also  makes  it  easy 
to  add  components  to  the  configuration,  to  move  two  bodies  relative  to  one  another, 
and  to  add  higher  resolution  grids  in  areas  of  interest. 

Blanking  or  I-blanking  as  it  is  sometimes  called  is  used  to  avoid  updating  the 
solution  inside  of  holes.  These  holes  come  about  in  several  ways.  The  most  obvious 
way  is  when  part  of  an  embedded  grid  lies  outside  of  the  domain  of  interest  (as  occurs 
when  a  cell  is  inside  of  a  solid  boundary).  This  case  is  not  encountered  in  this  study 
but  is  mentioned  for  completeness.  Holes  are  also  caused  when  two  or  more  grids 
cover  the  same  region  in  the  domain.  This  is  done  to  prevent  the  large  number  of 
interpolations  that  would  be  necessary  if  all  of  the  points  in  the  overlapping  region  of 
one  of  the  grids  were  to  be  set  by  the  other  overlapping  grid.  Not  only  could  a  large 
number  of  interpolations  be  expensive,  they  could  also  degrade  the  global  accuracy 
when  cells  of  the  overlapping  grids  are  different  sizes  [5].  In  this  case  the  grid  that 
provides  the  best  physical  solution  should  be  used  for  the  solution  procedure  in  the 
overlap  region.  This  can  be  difficult  to  determine  in  most  cases  and  is  usually  up  to 
the  engineer  to  decide  which  grid  will  yield  the  best  solution. 

Once  a  point  is  determined  to  be  in  a  hole  or  on  a  hole  boundary  it  becomes 
blanked.  The  residuals  calculated  for  these  blanked  points  are  set  to  zero  so  that 
no  correction  is  made  to  the  solution.  Since  no  corrections  are  made  to  the  solution 
for  blanked  points  boundary  conditions  need  to  be  applied  on  the  hole  boundaries. 
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This  is  accomplished  by  interpolating  the  dependent  variables  from  the  grid  that  is 
causing  the  blanked  points.  Grid  boundaries  that  lie  within  another  grid  must  also 
obtain  boundary  information  in  this  way.  This  requires  interpolation  stencils  to  be 
found  in  the  grid  causing  the  hole  for  these  boundary  points.  It  can  be  difficult  to 
find  valid  interpolation  stencils.  Valid  stencils  do  not  use  points  that  are  themselves 
set  by  interpolation.  In  the  case  where  no  valid  stencil  can  be  found  the  point  in 
question  is  referred  to  as  an  “orphan  point”  [121.  This  requirement  usually  leads  to 
an  overlap  of  two  to  three  cells  depending  on  the  interpolation  stencil.  Best  results 


tend  to  be  obtained  when  grid  cells  are  aoproximately  the  same  size  in  overlapping 
grids  [5,  12]. 

In  order  to  carry  out  the  intergrid  communication  a  method  of  finding  the  interpo¬ 
lation  stenciis  must  be  employed.  The  “stencil  jumping”  or  “stencil  walking”  [5.  12] 
technique  is  a  very  efficient  method  of  finding  the  nearest  neighboring  point.  Since 
the  FDS  algorithm  is  a  finite  volume  method  (ie.  the  dependent  variables  are  stored 
at  the  cell  centers)  the  cell  centered  meshes  must  be  constructed  in  order  to  search 
for  interpolation  stencils. 

This  stencil  jumping  technique  is  actually  an  inverse  transfinite  interpolation  prob¬ 
lem.  By  assuming  a  linear  variation  in  ^  and  t)  between  0  and  1  in  each  interpolation 
stencil  with  cell  centers  as  corners  Eq.  43  can  be  obtained. 
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and  is  the  lower  left  corner  point  of  the  cell  being  searched.  The  point 

A'p,  is  the  cell  center  from  another  grid  that  requires  an  interpolation  stencil.  The 
derivatives  in  the  matrix  in  Eq.  43  are  given  by 
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Eq.  43  is  in  the  form  of  a  Newton’s  method  and  will  converge  at  nearly  a  quadratic 
rate  as  AA'’  and  Al'*”’’  go  to  zero.  This  method  works  well  on  fairl}"  uniform  grids 
without  singularities. 

Once  the  nearest  neighbor  is  found  interpolation  can  be  applied  to  obtain  values 
at  the  boundary  points.  This  study  uses  bi-linear  interpolation  given  by 
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This  interpolation  procedure  in  general  does  not  maintain  conservation  principles  but 
it  is  commonly  used  in  this  fashion  with  satisfactory  results  [5].  The  use  of  conser¬ 
vative  techniques  (such  as  those  in  references  [6,  20]  )  would  require  a  substantial 
programming  effort  in  order  to  determine  areas  of  intersection  between  overlapping 
cells.  This  would  also  complicate  interpolation  stenciles  and  can  cause  stability  prob¬ 
lems. 

Note  that  in  the  equations  in  this  section  {i,j)  denotes  a  general  point  on  the  cell 


centered  grid. 


CHAPTER  5 

NONALIGNED  MULTIGRID  (NAM) 


The  idea  here  is  lo  treat  the  overlapping  grids  of  the  embedded  approach  as  in¬ 
dependent  levels  of  a  multigrid  FAS  procedure.  As  with  the  FAS  scheme,  defect 
corrections  and  dependent  variables  are  restricted  to  the  lower  levels  and  corrections 
are  prolongated  to  the  upper  levels.  Unlike  the  standard  FAS  where  one  level  directly 
communicates  with  only  the  next  coarser  level,  NAM  levels  may  directly  restrict  or 
prolongate  information  to  any  or  all  levels  depending  on  the  particular  system  of 
grids.  NAM  is  actually  a  generalized  local  refinement  method  [9,  11,  21]  in  the  sense 
that  the  various  levels  are  in  general  not  related  and  may  not  have  any  information 
restricted  to  portions  of  their  domain.  The  lack  of  relationship  between  levels  cov¬ 
ering  the  same  domain  is  similar  to  the  multigrid  methods  applied  to  unstructured 
grids  [20].  Unlike  embedded  grids  where  blanked  hole  points  are  eliminated  from  the 
solution,  NAM  uses  the  regions  of  overlapping  grids  to  generate  lower  level  corrections 
producing  the  faster  convergence  of  a  multilevel  scheme.  This  is  possible  because  the 
defect  correction  stops  grid  cells  of  different  sizes  from  degrading  the  global  solution. 
Blanking  is  still  appropriate  in  NAM  for  cases  where  the  cell  centers  of  one  grid  level 
are  located  inside  of  a  solid  surface  on  another  level  or  when  a  portion  of  one  of  the 
levels  is  outside  of  the  domain  of  interest.  Problems  with  blanked  points  in  the  lower 
levels  of  a  multigrid  procedure  are  discussed  in  references  [3]  and  [17]. 

The  communications  between  l^*veis  can  be  illustrated  by  applying  the  NAM  pro¬ 
cedure  to  the  system  of  embedded  grids  depicted  in  Figure  2.  Treating  this  system 
with  the  NAM  procedure  the  levels  are  as  indicated,  assuming  that  level  1  provides 
the  best  resolution  and  that  level  2  provides  better  resolution  than  level  3.  In  this 
case  Level  1  restricts  information  to  level  2  in  region  1  and  the  remainder  of  level 
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1  restricts  information  to  level  3.  Level  2  restricts  information  to  level  3.  The  re¬ 
stricted  information  in  region  1  contains  information  from  level  1  in  the  level  2  defect 
correction.  On  the  way  back  up  level  3  prolongates  the  lower  level  correction  to  all 
of  level  2  (including  region  1)  and  to  level  1  except  in  region  1.  Level  2  prolongates 
a  correction  to  level  1  in  region  1.  A  grid  schedule  analogous  to  the  FAS  V-Cycle 
diagram  is  shown  in  Figure  4  for  the  three  level  treatment  of  grids  just  described. 

Another  interesting  aspect  of  local  refinements  in  general  is  having  multiple  grids 
on  the  same  level  as  indicated  in  Figure  3  where  level  1  contains  two  grids  as  indicated. 
This  grid  is  treated  as  part  of  level  1  because  it  receives  no  information  from  a  higher 
level.  When  this  occurs  the  numerical  algorithm  is  applied  to  both  grids  independently 
and  information  is  restricted  to  levels  2  and  3  as  required.  In  general,  multiple  grids 
can  occur  on  every  level  and  when  coding  the  NAM  procedure  no  designation  needs 
to  be  made  between  levels  and  grids.  This  means  that  grid  2  could  be  treated  as  just 
another  level  without  any  prolongation  or  restriction  with  non-overlapping  levels. 

For  faster  convergence  of  this  problem  any  of  the  grids  may  be  coarsened  and 
added  as  additional  levels.  Care  should  be  used  when  coarsening  levels  other  than 
the  lowest  level.  For  example,  coarsening  level  1  and  inserting  ii  as  level  2  may  result 


Figure  2;  NAM  Three  Level  System 
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in  a  grid  with  worse  resolution  than  the  level  receiving  the  restricted  information. 
This  would  result  in  forcing  lower  resolution  onto  the  lower  levels.  The  lowest  level 
(typically  the  global  grid)  is  usually  the  safest  to  coarsen.  Coarsening  the  global  grid 
will  in  general  provide  faster  convergence  by  allowing  larger  time  steps  for  a  more 
rapid  development  of  the  flow  structure. 


Figure  3;  NAM  Three  Level  System  with  Multiple  Grids  per  Level 


Restriction  and  Prolongation  Operators 


As  in  Chapter  3  these  operators  define  the  interpolation  needed  to  move  quantities 
between  grids.  The  arbitrary  orientations  permissible  with  NAM  require  a  more  gen¬ 
eral  set  of  operators.  The  interpolation  stencils  are  found  using  the  “stencil  jumping” 
technique  discussed  in  Chapter  4.  Bi-linear  interpolation  is  used  for  both  restriction 
and  prolongation  operators  as  indicated  in  Eqs.  47,  48,  and  49.  The  subscripts  {i,j) 
denote  a  general  point  in  the  higher  level  and  {k,l)  denotes  a  general  point  in  the 
lower  level. 
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In  tne  previous  equcition  'l?  is  the  volume.  The  appearance  of  the  volume  is  due  to 
treating  the  residuals  as  point  quantities  during  the  interpolation  process. 
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The  bi-linear  interpolation  is  not  conservative,  but  is  selected  here  based  on  its  abil¬ 
ity  to  handle  communications  with  embedded  grids  [5j.  In  Chapter  6  FAS  produces 
the  same  results  as  the  single  grid  calculation  using  nonconservative  prolongation 
and  restriction  operators.  These  results  combined  with  successful  applications  on 
embedded  grids  would  tend  to  indicate  that  bi-linear  interpolation  is  worth  trying. 
Alternatives  include  the  conservative  techniques  used  in  unstructured  multigrid  '20' 
and  the  techniques  discussed  in  references  [6,  7]. 


V-Cycle  Procedure 


Level 


Figure  4:  NAM  Three  Level  Grid  Schedule 


The  cycling  procedure  and  grid  scheduling  is  similar  to  that  of  the  FAS  method 
and  FAS  is  actually  a  special  case  of  NAM.  The  NAM  grid  schedule  shown  in  Figure  4 
is  for  a  three  level  configuration  as  in  Figure  2.  Several  methods  of  scheduling  are 
possible  depending  on  the  system  of  grids  to  be  used.  The  solution  procedure  follows 
for  the  grid  schedule  in  Figure  4. 
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1.  Iterate  i/i  times  on  =  —At(L^{q^)) 
for  (  n  =  2;  n  <  3;  n  =  n+1  )  { 

2.  Calculate  the  residual  =  L''~^{q^~^) 
for  (  m  =  n+1;  m  <  3;  m  =  m+1  )  { 

3.  Restrict  solution  to  lower  level  q'^  = 

4.  Restrict  residual  minus  defect  correction  to  lower  level 
(recall  =  0) 

} 

5.  Calculate  the  defect  correction  t”  =  —  R” 

6.  Iterate  i/„  times  on  Aq"  =  — Af(I-”(g”)  —  t") 

} 

for  (  n  =  3;  n  >  2;  n  =  n-1  )  { 

7.  Calculate  the  lower  level  correction(LLC)  F"  =  Q"  — 
for  (  m  =  n-1;  m  >  1;  m  =  m-1  )  { 

8.  Prolongate  LLC  and  add  to  solution  on  finer  level  7”*  =  9™  +  /^V'" 

} 

9.  Iterate  i/„  times  on  Ag"“^  =  — 

} 


10.  If  solution  is  not  converged  go  to  1 


CHAPTER  6 
RESULTS 


The  two  dimensional  5°  ramp  near  a  flat  plate  shown  in  Figure  5  provides  the  test 
case  for  comparing  the  techniques  discussed  in  Chapters  2  through  5.  By  accurately 
predicting  the  flow  field  for  this  test  configuration  the  feasibility  of  the  nonaligned 
multigrid  technique  is  shown.  Solutions  generated  by  the  rest  of  the  methods  discussed 
are  to  show  the  relative  improvement  of  the  nonaligned  multigrid  and  to  provide  an 
understanding  of  how  the  various  components  of  the  nonaligned  multigrid  work.  The 
algorithms  are  coded  in  “C”  and  the  numerical  solutions  are  generated  using  double 
precision  numbers  on  an  Iris  4D/320VGX  workstation.  No  attempt  was  made  to 
parallelize  on  this  two  processor  machine. 

Theoretical  Solution 


The  theoretical  solution  to  this  test  case  is  obtained  from  the  oblique  shock  charts 
and  the  normal  shock  tables  for  isentropic  flow  [19].  The  weak  shock  solution  is  given 
in  Table  1  and  the  shock  angles  are  depicted  in  Figure  5. 


Table  1:  Theoretical  Solution 
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Region 

2.00 

P 

1.00 

pu 

2.00 

pv 

0.00 

e 

3.79 

P 

0.71 

I  Region  II 
1.83 
1.20 
2.28 
0.20 
4.49 
0.93 


Region  III 
1.60 
1.49 
2.58 
0.00 
5.37 
1.25 
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Figure  5;  Problem  Geometry 


Single  Grid  Calculations 

A  single  grid  solution  is  obtained  by  applying  the  flux  difference  splitting  (FDS) 
algorithm  to  the  40  by  32  grid  in  Figure  6.  Density  contours  from  the  solution  are 
plotted  in  Figure  7  after  converging  to  machine  zero.  This  solution  does  show  the 
general  features  of  the  shock  structure;  however,  the  first  order  FDS  and  the  fairly 
coarse  grid  combine  to  smear  the  shock  waves  and  the  zero  pressure  gradient  boundary 
conditions  cause  the  shocks  to  curve  near  the  solid  surfaces.  In  order  to  gain  better 
accuracy  grids  with  dimensions  80  by  64  and  160  by  128  are  also  used  to  generate 
solutions.  The  80  by  64  grid  plotted  in  Figure  8  is  obtained  by  quartering  the  grid 
cells  in  the  40  by  32  grid.  The  160  by  128  grid  is  too  dense  to  be  plotted  clearly  and 
is  obtained  by  quartering  the  80  by  64  grid  cells.  Density  contours  of  the  resulting 
solutions  are  shown  in  Figures  9  and  10.  The  contour  plots  show  that  as  the  number 
of  grid  points  is  increased  the  dissipation  becomes  less  and  the  curving  near  the  solid 
boundaries  decreases  causing  the  solution  to  come  closer  to  the  theoretical  values  in 
Table  1. 

The  next  area  of  interest  is  to  determine  the  efficiency  of  the  solution  procedure 
and  the  overall  cost.  In  order  to  do  this  the  following  measure  of  the  change  in  the 
solution  is  generated  for  each  iteration. 


R  =  logjo 


jmax  imai  T  4 

E  E  lE 

j=i  1=1  U=i 

By  looking  at  the  history  of  this  residual,  plotted  in  Figure  11.  and  the  CPU  times 
in  Table  2  the  cost  of  using  an  explicit  scheme  becomes  apparent.  The  explicit  FDS 
scheme  converges  in  a  fairly  slow  manner  due  to  the  stability  condition  indicated  in 
Eq.  24.  Since  this  is  a  steady  state  problem  local  tirre  stepping  is  applied  with  a 
Courant  number  of  0.96  in  order  to  provide  the  largest  time  step  in  each  grid  cell. 
With  the  addition  of  grid  points  the  already  slow  convergence  rate  becomes  even 
slower  due  to  the  stability  condition  requiring  smaller  time  steps  as  the  grid  spacing 
decreases.  This  addition  of  grid  points  not  only  decreases  the  convergence  rates,  but 
it  also  increases  the  CPU  time  as  indicated  in  Table  2.  Although  these  differences 
in  CPU  time  and  convergence  rates  may  not  seerri  that  important  for  a  problem  of 
this  size,  these  factors  become  very  limiting  for  more  complicated  three  dimensional 
problems. 
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Multigrid  Calculations 


The  full  approximation  scheme  (FAS)  discussed  in  Chapter  3  is  applied  to  the 
g'ids  of  the  previous  section  in  order  to  obtain  faster  convergence. 

Grid  Coarsening 


The  first  step  is  to  coarsen  the  grids  as  indicated  in  Figure  12  for  the  40  by  32 
grid.  The  coarsening  used  here  simply  removes  every  other  grid  point  on  the  current 
level  in  each  curvilinear  coordinate  direction  to  produce  the  next  coarser  grid  which 
will  contain  1/4  the  number  of  points  of  the  current  level  [9,  10].  When  doing  this 
it  is  important  to  preserve  the  geometry  on  the  lower  levels.  This  is  accomplished 
here  by  placing  a  sufficient  number  of  grid  points  in  front  of  the  ramp  such  that  the 
leading  edge  point  will  not  be  eliminated  by  the  coarsening  procedure  for  the  number 
of  levels  that  are  used  in  these  calculations. 

Work  Units 


Before  making  any  calculations  the  multigrid  work  unit  [2,  9,  10]  needs  to  be 
defined.  The  work  unit  is  the  ratio  of  the  number  of  grid  points  on  the  current  level 
to  the  number  of  grid  points  on  the  fine  level.  This  means  the  work  units  are  ]  on  level 
h,  1/4  on  level  2h,  and  1/16  on  level  4h  for  the  grid  coarsening  discussed  above.  With 
this  definition  the  work  required  for  the  restriction  and  prolongation  operations  is 
neglected.  This  is  the  standard  multigrid  work  unit  definition  which  seems  reasonable 
since  the  cost  of  interpolation  between  levels  is  insignificant  compared  to  the  cost  of 
the  FDS  solution  procedure.  With  this  definition  the  work  unit  provides  a  means  of 
comparison  between  the  single  grid  and  multigrid  convergence  histories. 
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FAS  Results 


Applying  the  FAS  procedure  to  the  grids  (40  by  32,  80  by  64,  and  160  by  128)  dis¬ 
cussed  in  the  previous  section  yields  the  same  solution  (to  double  precision  accuracy) 
as  those  in  Figures  7,  9,  and  10,  respectively.  This  result  would  seem  to  indicate  that 
the  nonconservative  restriction  and  prolongation  operators  defined  in  Chapter  2  are 
accurate  enough  to  keep  from  introducing  significant  error  into  the  solution. 

The  residual  histories  for  the  two  and  three  level  calculations  on  the  40  by  32 
grid  are  plotted  in  Figure  13  along  with  the  single  grid  residual.  Several  calculations 
were  made  varying  the  iteration  count  on  each  level  in  order  to  find  the  maximum 
gain  in  convergence.  It  was  fairly  easy  to  find  a  good  two  level  calculation,  but 
when  an  additional  level  is  added  the  possible  number  of  variations  becomes  large. 
The  iteration  counts  used  for  Figure  13  are  given  in  Table  3.  Comparing  the  two 
and  three  level  residuals  indicates  that  the  three  level  calculation  could  not  achieve 
an  increase  in  performance  over  the  two  level  calculation.  Also  note  that  the  best 
convergence  on  the  three  level  calculation  was  achieved  by  running  one  iteration  per 
multigrid  cycle  on  the  coarsest  level.  The  combination  of  these  two  would  seem  to 
indicate  that  the  4h  level  is  too  coarse  to  handle  the  shock  wave.  The  wider  cells  on 
the  4h  level,  depicted  as  the  bottom  grid  in  Figure  12,  are  more  likel}’  to  cause  the 
restriction  and  prolongation  operators  to  straddle  the  shock  wave  resulting  in  faulty 
communication  between  levels  [9].  A  conservative  set  of  operators  should  help  to 
overcome  this  problem. 

Similar  results  are  obtained  by  looking  at  the  convergence  histories  for  the  80  by  64 
grid  in  Figure  8.  Figure  14  plots  the  residual  histories  for  the  single  grid  calculation 
and  for  the  two  and  three  level  FAS  calculations.  In  this  case  the  convergence  is 
slightly  better  for  the  three  level  calculation.  The  residual  history  indicates  that  the 
third  level  helps  to  accelerate  convergence  early  in  the  calculations,  bui  the  gain  is 
lost  as  the  solution  becomes  more  accurate.  No  four  level  calculations  are  attempted 
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for  this  case  since  the  three  level  results  on  the  40  by  32  grid  did  not  provide  any 
improvement.  The  iteration  counts  used  in  these  calculations  are  given  in  Table  3. 

Figure  15  shows  the  residual  histories  for  the  single  grid  and  the  two,  three,  and 
four  level  calculations  on  the  160  by  128  grid.  In  this  case  the  four  level  V-Cycle 
provides  the  best  results  as  expected. 

Several  observations  can  be  made  from  the  convergence  histories.  First,  the  best 
results  occur  when  smoothing  iterations  are  made  (i/  7=  0)  when  prolongating  the 
coarse  grid  correction.  Secondly,  the  40  by  32  grid  results  show  that  there  is  a  practical 
limit  as  to  how  course  the  grids  can  be  and  still  help  to  reduce  the  number  of  work 
units  for  the  hyperbolic  problem  when  nonconservative  interpolation  operators  are 
used.  Additionally,  a  method  of  automatically  determining  the  iteration  count  on  the 
various  levels  [9,  17]  needs  to  be  implemented  based  on  the  number  of  runs  made  to 
find  the  iteration  counts  for  the  calculations  presented  here. 

A  direct  consequence  of  the  faster  convergence  rate  provided  by  multigrid  is  the 
reduced  CPU  time.  As  noted  previously  the  restriction  and  prolongation  operations 
were  neglected  in  the  work  unit  calculations  but  are  included  in  the  CPU  time.  This 
explains  why  the  savings  in  CPU  time  indicated  in  Tables  3,  4,  and  5  may  not  be  as 
large  as  the  savings  indicated  by  the  work  units. 

Table  3:  FAS  Calculation  Statistics  on  the  40  X  32  Grid 

Levels  i/h  U2h  ^2h  i^h  Work  Units  Multigrid  Cycles  CPU  Time 

1  1  —  _  _  _  550.000  —  131.6  sec 

2  2  6  -  -  2  396.000  61  94.8  sec 

3  3  5  1  4  2  411.000  48  96.3  sec 
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Table  4:  FAS  Calculation  Statistics  on  the  80  X  64  Grid 


Levels  Uh  i>2h 

1  1  -  - 

2  2  4  - 

3  4  4  14 


^2h  Work  Units  Multigrid  Cycles  CPU  Time 

—  —  855.000  —  900.2  sec 

-  2  552.000  92  547.9  sec 

4  4  545.625  45  547.0  sec 


Table  5:  FAS  Calculation  Statistics  on  the  160  X  128  Grid 


Levels 

^2h 

i'2h 

i>h 

Work  Units 

Multigrid  Cycles 

CPU  Time 

1 

1 

— 

— 

— 

— 

— 

— 

1412.00 

— 

5863.0  sec 

2 

4 

8 

— 

— 

— 

— 

5 

804.00 

61 

3345.3  sec 

3 

10 

8 

18 

— 

— 

8 

10 

791.25 

30 

3316.9  sec 

4 

12 

9 

8 

14 

8 

9 

12 

775.78 

25 

3034.7  sec 

Figure  12:  40  X  32  Grid  and  Two  Coarser  Levels 


Log(L2  -  NORM)  Log(L2 


Figure  14:  FAS  Residual  Histories  on  80  X  64  Grid 
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Figure  15:  FAS  Residual  Histories  on  160  X  128  Grid 


Embedded  Grid  Calculations 


The  embedded  grid  approach  is  used  here  to  provide  an  improved  solution  by 
inserting  local  grids  over  the  regions  where  the  shock  waves  occur.  Several  variations 
are  attempted  in  order  to  predict  the  theoretical  solution  as  accurately  as  possible. 

5°  Ramp  Calculations 


In  order  to  provide  a  better  understanding  the  embedded  approach  is  first  imple¬ 
mented  on  a  5°  ramp.  This  eliminates  the  reflected  shock  and  makes  the  results  easier 
to  interpret.  A  16  by  16  local  grid  is  embedded  into  a  40  by  32  grid  as  indicated  in 
Figure  16.  The  cell  centered  grids  are  plotted  in  Figure  17  to  show  the  hole  caused  by 
the  local  grid.  Since  all  boundary  points  are  blanked  the  overlapping  regions  appear 
to  be  much  narrower  than  they  actually  are. 
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Aligning  the  local  grid  lines  parallel  to  the  shock  wave  allows  the  FDS  to  produce 
a  much  more  accurate  solution  than  is  possible  on  the  global  grid.  The  reason  for 
the  improved  accuracy  is  due  to  the  assumption  in  FDS  that  the  characteristic  waves 
propagate  in  a  direction  perpendicular  to  the  grid  interfaces  [32]. 

The  density  contours  from  this  two  grid  solution  are  plotted  in  Figure  18.  This 
solution  is  much  better  than  could  be  obtained  on  the  global  grid  alone  and  is  virtually 
the  same  as  the  theoretical  solution  (  regions  I  and  II  in  Table  1  and  in  Figure  5). 

Figure  19  shows  the  residual  histories  for  the  global  grid  alone  and  for  the  two 
grid  case.  The  residual  for  the  embedded  grid  is  obtained  by  using  Eq.  50  on  each 
grid.  For  blanked  points  no  contribution  is  made  to  the  residual.  This  plot  shows 
that  the  cost  of  the  embedded  grid  calculation  (340.4  seconds)  is  greater  than  for 
the  single  grid  calculation  (96.2  seconds);  however,  the  solution  on  the  embedded 
system  is  so  superior  to  that  of  the  single  grid  that  this  cost  comparison  is  virtuali}' 
meaningless.  For  the  single  grid  calculation  the  dissipation  of  the  first  order  FDS 
causes  the  width  of  the  predicted  shock  to  be  slightly  wider  than  the  width  of  the 
local  grid  in  Figure  16  at  the  back  plane.  The  increased  accuracy  of  this  embedded 
solution  becomes  apparent  when  comparing  the  width  of  the  smeared  shock  to  that 
of  the  shock  predicted  in  Figure  18. 


Figure  16:  Embedded  Grids  for  the 
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Figure  18:  Density  Contours  for  the  5°  Ramp  Embedded  Solution 


Figure  19:  Residual  Histories  for  the  5°  Ramp  Embedded  Calculation 
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5°  Ramp  Near  a  Flat  Plate  Calculations 


Two  grid  case 

Figure  20  shows  the  embedded  grid  configuration  with  the  lines  of  the  local  grid 
again  parallel  to  the  shock.  The  hole  created  in  the  global  grid  can  be  seen  in 
Figure  21. 

Density  contours  of  the  solution  obtained  on  this  grid  system  are  shown  in  Fig¬ 
ure  22.  Notice  how  much  more  accurately  the  initial  shock  is  captured  than  in  the 
single  grid  solution  (compare  to  Figure  7).  The  reflected  shock  is  not  resolved  due  to 
the  alignment  of  the  local  grid  lines. 

Comparing  the  residual  histories  plotted  m  Figure  23  with  those  in  Figure  11 
for  the  40  by  32  single  grid  calculation  shows  that  the  number  of  iterations  is  more 
than  double  that  of  the  single  grid  at  machine  zero.  A  large  part  of  this  is  due 
to  the  introduction  of  the  artificial  boundaries  into  the  domain  (the  front  and  back 
boundaries  of  the  local  grid).  The  explicit  communication  between  grids  at  these 
boundaries  hinders  the  development  of  the  solution. 

The  solution  on  this  system  of  grids  requires  425.3  CPU  seconds  to  converge  to 
machine  zero.  This  is  more  than  three  times  the  cost  of  the  single  grid  calculation. 
This  increase  reflects  the  cost  of  making  additional  calculations  on  the  local  grid,  the 
cost  of  interpolations  between  grids,  and  the  cost  of  running  more  iterations  due  to 
the  slower  convergence. 
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Figure  22:  Two  Embedded  Grid  Density  Contours 


Figure  23:  Two  Embedded  Grid  Residual  Histories 
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Three  grid  case 


Based  on  the  results  of  the  two  grid  calculations  a  third  grid  that  is  aligned  with 
the  reflected  shock  is  inserted  as  indicated  in  Figure  24.  The  plot  of  the  cell  centered 
grids  in  Figure  25  shows  the  holes  for  this  grid  system.  In  this  case  the  local  grid 
from  the  previous  section  (referred  to  as  grid  1)  punches  holes  in  both  the  new  local 
grid  (grid  2)  and  in  the  global  grid  while  grid  2  only  punches  a  hole  in  the  global  grid. 

The  solution  obtained  with  this  grid  arrangement  is  shown  in  Figure  26.  The  solu¬ 
tion  remains  good  for  the  initial  shock  but  the  reflected  shock  is  still  not  sufficiently 
resolved.  Since  the  reflected  shock  is  not  resolved  on  grid  1  it  cannot  feed  good  values 
to  grid  2  hence  grid  2  is  incapable  of  capturing  the  reflected  shock  accurately. 

Residual  histories  for  this  calculation  are  in  Figure  27.  Once  again  the  cost  is 
increased  by  the  addition  of  the  embedded  grid.  This  grid  system  requires  514.9 
CPU  seconds  to  reach  machine  zero. 


Figure  24:  Three  Embedded  Grids  for  the  5°  Ramp  Near  a  Flat  Plate 


Figure  25:  Three  Grid  Blanking  for  the  5°  Ramp  Near  a  Flat  Plate 
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Figure  26:  Three  Embedded  Grid  Density  Contours 
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Four  grid  case 


The  addition  of  a  third  local  grid  (grid  3)  is  made  in  hopes  of  accurately  repre¬ 
senting  both  shocks  near  the  reflection  point.  Grid  3  is  a  semicircle  with  its  center  at 
the  shock  reflection  point  as  indicated  in  Figure  28.  The  radial  lines  in  this  grid  are 
arranged  such  that  one  grid  line  will  coincide  with  the  initial  shock  and  another  line 
with  the  reflected  shock.  In  using  this  grid  a  singularity  is  introduced  at  the  shock 
reflection  point.  When  the  FDS  scheme  is  applied  to  this  grid  the  fluxes  at  this  point 
are  set  to  zero  to  prevent  flow  through  the  solid  wall.  This  new  grid  punches  holes  in 
the  three  previous  grids  as  shown  in  Figure  29. 

Density  contours  of  the  solution  obtained  for  this  system  are  plotted  in  Figure  30. 
This  new  solutiori  is  very  accurate  and  is  almost  indistinguishable  from  the  theoretical 
solution. 

The  residual  history  for  this  grid  system  is  plotted  in  Figure  31  and  the  CPU 
time  to  reach  this  point  is  865.0  seconds.  Although  this  embedded  system  requires 
many  times  the  amount  of  work  needed  for  the  single  grid,  the  solution  obtained  here 
is  far  superior  to  that  in  Figure  7.  Actually,  this  embedded  solution  is  better  than 
the  solution  generated  on  the  160  by  128  single  grid  and  costs  less  to  produce.  The 
iterations  in  Figure  11  must  be  multiplied  by  16  to  compare  with  the  embedded  grid 
work  units  (this  is  due  to  one  embedded  work  unit  being  equal  to  one  iteration  on 
the  40  by  32  grid). 


Figure  28:  Four  Embedded  Grids  for  the  5°  Ramp  Near  a  Flat  Plate 


Figure  29:  Four  Grid  Blanking  for  the  5°  Ramp  Near  a  Flat  Plate 
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Figure  30:  Four  Embedded  Grid  Density  Contours 
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Figure  31:  Four  Embedded  Grid  Residual  Histories 
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Nonaligned  Multigrid  Calculations 


The  goal  of  this  section  is  to  combine  the  improved  solution  obtained  in  the  em¬ 
bedded  approach  with  the  improved  convergence  of  the  multigrid  in  order  to  show 
the  feasibility  of  the  nonaligned  multigrid  approach. 

Restriction  and  Prolongation  Operators 


Prior  to  applying  the  nonaligned  multigrid  the  restriction  and  prolongation  oper¬ 
ators  need  to  be  changed  as  indicated  in  Chapter  5.  In  order  to  determine  how  well 
these  new  operators  perform  one  of  the  two  level  FAS  calculations  was  repeated.  The 
same  solution  was  achieved  requiring  the  same  number  of  work  units  for  both  sets  of 
operators.  Since  the  same  solution  was  obtainable  in  the  same  number  of  work  units 
it  appears  that  conservation  is  no  more  of  a  problem  for  the  new  operators  than  for 
the  old  operators;  however,  this  does  not  indicate  what  will  happen  when  the  grids 
are  no  longer  aligned.  The  new  operators  required  approximately  two  seconds  more 
than  the  old  operators  for  the  test  case  (the  two  level  2-4-2  V-Cycle  on  the  40  by  32 
grid).  This  slight  increase  in  cost  is  attributed  to  the  time  used  to  apply  the  more 
complicated  interpolations  and  the  time  required  for  the  “stencil  jumping”  procedure. 

5°  Ramp  Calculations 


Here  the  NAM  procedure  is  applied  to  the  embedded  grid  system  of  Figure  16  to 
provide  accelerated  convergence.  Designating  the  local  grid  as  level  1  and  the  global 
grid  as  level  2  the  solution  in  Figure  32  is  obtained. 

Several  calculations  were  made  varying  the  iteration  counts  on  the  two  levels  in 
an  attempt  to  achieve  optimal  convergence.  The  resulting  number  of  work  units 
required  to  reach  machine  zero  are  listed  in  Table  6  for  five  of  these  runs.  This 
shows  the  difficulty  of  finding  the  best  iteration  count  per  level  to  provide  maximum 
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convergence.  This  problem  was  also  discussed  in  the  FAS  results.  Since  the  severity 
of  this  problem  increases  with  the  addition  of  more  levels  an  automated  switching 
scheme  is  needed. 

Eq.  50  is  applied  to  the  highest  level  at  each  point  in  the  domain  to  obtain  the 
residuals.  This  is  accomplished  by  removing  any  point  from  the  residual  calculation 
that  receives  information  from  a  higher  level.  The  residual  histories  for  the  best  case 
from  Table  6,  for  the  embedded  grid  calculations,  and  for  the  single  grid  calculation 
are  plotted  in  Figure  33.  This  shows  that  the  NAM  calculations  require  approximately 
one  third  the  number  of  work  units  required  by  the  embedded  grid  procedure.  For  the 
run  plotted,  NAM  requires  137.6  CPU  seconds  compared  to  the  340.4  CPU  seconds 
required  by  the  embedded  approach.  This  is  a  59%  reduction  in  CPU  time. 

Three  attempts  are  made  to  further  increase  the  convergence  rates  for  this  set  of 
grids.  The  first  attempt  is  to  coarsen  level  1  to  produce  a  new  level  2  and  moving  the 
global  grid  to  level  3.  The  resulting  three  level  NAM  procedure  requires  524.8  work 
units  and  133.6  CPU  seconds  to  reach  machine  zero  for  the  best  iteration  count  per 
level  (i/i  =  2,1/2  =  0,1/3  =  3,1/2  =  l,i/i  =  2).  This  is  not  a  significant  improvement 
over  the  two  level  NAM  results.  Level  2  in  this  case  acts  as  a  filter  to  smooth  the 
solution  and  residual  being  passed  from  level  1  to  level  3.  The  same  solution  is 
obtained  as  in  the  two  level  calculation. 

A  second  three  level  NAM  procedure  is  defined  by  making  the  local  grid  level  1, 
the  global  grid  level  2,  and  coarsening  level  2  to  make  level  3.  This  system  takes 
only  469.7  work  units  requiring  114.3  CPU  seconds  to  converge  to  machine  zero.  The 
iteration  count  used  here  is  vj  =  3,  t/2  =  2,1/3  =  2, 1/2  =  1,  and  i/i  =  4.  This  system 
requires  only  69.7  more  work  units  (5.2  seconds)  than  the  single  grid  calculation  while 
producing  a  far  superior  solution.  This  is  a  69%  reduction  in  the  CPU  time  required 
by  the  embedded  grid  approach. 
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Figure  32:  Density  Contours  for  the  5°  Ramp  NAM  Solution 


The  final  attempt  to  converge  faster  than  the  single  grid  is  made  by  generating 
a  four  level  NAM  system  by  coarsening  level  3  in  the  first  three  level  procedure 
discussed.  The  best  case  obtainable  here  was  474.95  work  units.  This  may  be  caused 
by  the  speed  at  which  the  solution  is  being  generated  on  the  coarse  grid  being  hindered 
by  the  two  sets  of  artificial  boundaries  (the  left  and  right  sides  of  levels  1  and  2). 

These  results  seem  to  indicate  that  the  most  significant  gain  in  convergence  will 
come  from  simply  implementing  NAM  on  the  existing  embedded  grid  system.  It  also 
appears  to  be  better  to  coarsen  the  global  grid  to  provide  additional  convergence. 


Table  6:  NAM  Work  Units 

i/j  1/2  i'l  Work  Units 

1  3  2  563.4 

3  2  2  588.4 

3  3  3  558.2 

342  571.0 

2  3  3  549.6 


Log  (L2 
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5°  Ramp  Near  a  Flat  Plate  Calculations 


Two  grid  case 

This  calculation  uses  the  grids  of  Figure  20  with  the  local  grid  considered  to  be 
level  1  and  the  global  grid  to  be  level  2.  Density  contours  of  the  resulting  solution 
are  plotted  in  Figure  34.  This  shows  the  same  inability  to  capture  the  reflected  shock 
as  the  embedded  solution  (Figure  22). 

The  residual  plot  in  Figure  35  compares  the  residuals  from  the  two  grid  embedded 
calculations  with  the  residuals  of  the  current  calculation.  This  plot  shows  that  NAM 
requires  under  half  the  number  of  work  units  needed  for  the  embedded  calculation. 
This  could  be  improved  by  coarsening  the  global  grid  but  it  is  not  attempted  here 
since  the  solution  is  not  that  good.  The  CPU  time  for  the  NAM  is  174.8  seconds 
while  the  embedded  approach  requires  425.3  seconds.  The  iteration  count  for  this 
case  was  i/j  =  l,U2  =  2,  and  u\  =  1. 
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Figure  34:  NAM  2  Level  Density  Contours 
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Three  grid  case 


In  this  case  the  grids  of  Figure  24  are  used.  Level  1  contains  the  local  grid  covering 
the  initial  shock,  level  2  is  the  grid  covering  the  reflected  shock,  and  level  3  is  the 
global  grid.  Density  contours  of  the  resulting  solution  are  plotted  in  Figure  36.  The 
solution  is  similar  to  the  embedded  results  (Figure  26),  both  are  having  difficulty 
capturing  the  reflected  shock  accuratehc  Once  again  this  is  due  to  the  inability  of 
the  grids  to  accurately  model  the  physics  near  the  top  wall. 

The  residuals  for  this  calculation  and  the  three  grid  embedded  solution  are  plotted 
in  Figure  37.  Fhis  case  also  requires  approximately  half  the  number  of  work  units 
required  for  the  embedded  procedure.  NAM  requires  259.7  seconds  compared  to  514.9 
for  the  embedded  calculation. 
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Figure  36:  NAM  3  Level  Density  Contours 
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Four  grid  case 


For  this  case  the  NAM  is  applied  to  the  four  grid  system  in  Figure  28  treating 
the  circular  grid  as  level  1,  the  grid  aligned  with  the  initial  shock  as  level  2,  the  grid 
aligned  with  the  reflected  shock  as  level  3  and  the  global  grid  is  level  4.  The  four 
level  grid  scheduling  diagram  is  shown  in  Figure  38.  The  circled  numbers  indicate 
the  number  of  iterations  performed  on  each  level  for  a  single  NAM  cycle.  Density 
contours  of  the  resulting  solution  are  plotted  in  Figure  39.  These  are  similar  to  that 
obtained  by  the  embedded  system. 

Residual  histories  for  this  case  are  plotted  in  Figure  40.  Comparing  these  to  the 
embedded  grid  residuals  in  Figure  31  shows  faster  convergence  is  achieved  by  the 
NAM.  CPU  time  required  for  NAM  is  712.5  seconds  compared  to  the  865.0  seconds 
required  by  the  embedded  grids. 

To  further  increase  the  convergence  the  global  grid  in  this  four  level  system  is 
coarsened  to  produce  a  five  level  system.  The  grid  schedule  is  as  shown  in  Figure  41 
with  the  iteration  count  as  indicated  in  the  circles.  The  same  solution  is  obtained  in 
only  623.8  CPU  seconds.  The  residual  histories  are  as  indicated  in  Figure  42.  This 
is  27.9%  faster  than  the  equivalent  embedded  system. 
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Figure  38;  NAM  4  Level  Grid  Schedule 
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Figure  39:  NAM  4  Level  Density  Contours 
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Figure  40:  NAM  4  Level  Residual  Histories 
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CHAPTER  7 
CONCLUSIONS 


The  explicit  first  order  flux  difference  splitting  (FDS)  was  used  to  solve  the  equa¬ 
tions  governing  inviscid  fluid  flow  on  a  single  grid  for  the  5°  ramp  and  the  5°  ramp 
near  a  flat  plate.  Calculations  were  made  for  the  Mach  2  case  using  local  time  step¬ 
ping  at  a  Courant  number  of  0.96.  A  grid  refinement  study  was  also  performed  on 
the  5°  ramp  near  a  flat  plate  to  show  grid  dependency  of  the  solution  and  to  show 
the  increasing  >_Oot  as  the  grid  becomes  finer. 

The  multigrid  full  approximation  scheme  (FAS)  was  applied  to  this  nonlinear  prob¬ 
lem  to  provide  increased  convergence  rates  and  reduced  central  processing  unit  (CPU) 
time  requirements.  The  solution  was  obtained  independent  of  the  number  of  levels 
and  independent  of  the  number  of  iterations  performed  on  those  levels.  The  conver¬ 
gence  rate  was  determined  to  be  very  sensitive  to  the  number  of  iterations  performed 
on  the  various  levels. 

Several  systems  of  embedded  grids  were  implemented  to  provide  increased  accuracy 
near  shock  waves.  Embedded  grids  were  aligned  with  the  shocks  to  take  advantage  of 
the  excellent  shock  capturing  capability  of  the  FDS  scheme.  Although  this  method 
increased  accuracy  it  degraded  the  convergence  rate  and  increased  the  required  CPU 
time. 

The  nonaligned  multigrid  was  introduced  to  provide  increased  convergence  for  the 
embedded  grid  systems  by  treating  the  individual  grids  as  levels  in  the  multigrid  so¬ 
lution  procedure.  This  technique  was  able  to  converge  69%  faster  than  the  embedded 
grid  procedure  for  the  5°  ramp,  and  28%  faster  for  a  more  complex  reflected  shock 
case. 
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A  great  deal  of  work  is  still  required  to  make  the  method  into  a  “black  box”  en¬ 
gineering  tool.  Some  of  the  topics  include  automatic  determination  of  the  num.ber 
of  iterations  to  perform  on  each  level,  automated  grid  scheduling  for  optimal  perfor¬ 
mance,  and  determining  the  allowable  change  in  grid  spacing  between  levels.  Another 
very  important  topic  is  developing  a  conservative  procedure  for  communication  be¬ 
tween  levels. 

The  current  follow-on  plan  is  to  implement  the  NAM  procedure  in  a  viscous  three 
dimensional  higher  order  implicit  computational  fluid  dynamics  (CFD)  code  currently 
under  development. 


APPENDIX 

STABILITY  ANALYSIS 
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The  first  order  linear  wave  equation  has  the  same  form  as  the  Euler  equations 


written  in  matrix  form  (Eq.  1,  Chapter  2). 

^  M  4.  ^  =  0 

dt  d(  drj 


(51) 


In  this  case  q  is  the  dependent  variable  and  the  fluxes  are  =  aq  and  =  bq  where 
a  and  b  are  constants. 

Applying  the  flux  difference  splitting  (FDS)  described  in  Chapter  3  produces  the 
following  discretization,  which  is  simply  the  scalar  equivalent  of  Eq.  19. 


-  < 


At 


(52) 


Using  Eqs.  20,  22,  and  noting  that  the  left  and  right  eigenvectors  are  simply  reciprocals 
for  the  scalar  case,  the  numerical  fluxes  are  reduced  to 


if()7+i/2.j  —  (9t+i.j  9i,j) 

(7r,)M+l/2  =  + 


(53a) 

(53b) 


The  actual  fluxes  /*,  evaluated  at  the  discrete  points  are  also  represented  in  terms  of 
the  dependent  variable.  The  values  for  the  non-positive  eigenvalues  depend  on  a  and 
b  as  indicated  below.  The  case  of  a  or  fe  equal  to  zero  is  not  considered  since  it  is 
trivial  for  the  two  dimensional  case. 


Case  1  for  a  >  0  and  fe  >  0  — >•  =  0,  A~  =  0 

Case  2  for  a  <  0  and  6  >  0  — >  A^  =  a,  A“  =  0 

Case  3  for  a  >  0  and  6  <  0  — >  Aj  =  0,  A“  =  6 

Case  4  for  o  <  0  and  b  <  0  Xj  =  a,  X~  =  b 

Substituting  the  numerical  fluxes  into  Eq.  52  and  rearranging  yields 


(7"+i  =  <7" 


l<'(C  -  *-.)  +  +  *-.) 
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(54) 
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The  next  step  is  to  apply  the  von  Neumann  (or  Fourier)  stability  analysis  [1].  This 


is  done  by  substituting 


Where  i  =  ^  —1  into  Eq.  54.  Making  the  substitution,  using  trigonometric  relations 
and  grouping  terms  in  real  and  imaginary  parts  yields 


C  =  l  — 2At  (a  —  2Aj  )  sin'^(/S/2)  +  (6  —  2A^  )  sin^(Q/2)  —  zAt  [a  sin/3  +  6sin  a]  .  (56) 


Once  the  equation  is  in  this  form  require  |(^|  <  1  to  obtain  the  stability  condition 
on  the  time  step  size.  Applying  this  restriction  and  solving  for  At. 


At  < 


4  f(a  —  2A^  )sin^(/3/2)  +  (6  —  2A^  )sin^(Q/2)] 


4  (a  —  2A^  ) sin^(/S/2) -h  (6  —  2A“ ) sin^(Q:/2)l  +  [asin/3  +  isina 


The  four  possible  variations  on  the  eigenvalues  listed  previously  are  now  applied 
to  this  equation  and  the  most  restrictive  values  for  q  and  0  are  sought. 

Case  1;  Eq.  57  becomes 

4  (asin2(^/2)  +  fesin2(a/2)) 

4  \  as\T0{0 j2)  +  6sin^(a/2)j  +  (asin/3  +  6sin  a)^ 

The  minimum  value  of  the  right  hand  side  of  this  equation  occurs  at  a  = 

/3  =  TT  causing  the  stability  condition  to  become 


At  < 


a  +  6 


Case  2:  Eq.  57  becomes 

4  f— asin^(/3/2)  +  i)sin^(Q/2)) 

At  <  - - — — - 2 - - .  (60) 

4  a sin^(/3/2)  +  6sin^(Q/2)j  +(  asin/3  +  isina)^ 

The  minimum  value  of  the  right  hand  side  of  this  equation  occurs  at  a  = 

/3  =  TT  causing  the  stability  condition  to  become 


At  < 


b  —  a 


(61) 
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Case  3;  Eq.  57  becomes 


At  < 


4  ^asin^(/3/2)  —  6sin^(a/2)j 
4  ^asin^(/3/2)  —  fesin^(a/2)^  +  (asin/3  +  isina)^ 


The  minimum  value  of  the  right  hand  side  of  this  equation  occurs  at  a 
(3  —  TT  causing  the  stability  condition  to  become 


At  < 


a  —  b 


Case  4:  Eq.  57  becomes 


At  < 


4  (^— asin^(/3/2)  —  isin^(a/2)) 

4  (^—asm^{/3/2)  —  bsm^{a/2)J  -I-  (a sin/3  +  /)sin a)^ 


The  minimum  value  of  the  right  hand  side  of  this  equation  occurs  at  q  = 
/3  =  TT  causing  the  stability  condition  to  become 


At  < 


—  a  —  b 


The  stability  condition  in  general  for  FDS  algorithm  applied  to  the  two  dimensional 
wave  equation  can  be  written  as 


(66) 

l“l  +  161 

It  should  be  noted  that  this  analysis  assumes  a  linear  difference  equation  with 
periodic  boundary  conditions.  The  standard  practice  is  to  ignore  these  assumptions 
which  usually  provides  good  results.  To  apply  this  condition  to  the  Euler  equations 
the  maximum  eigenvalues  (Eq.  14,  Chapter  2)  of  the  flux  Jacobian  matrix  are  sub¬ 
stituted  into  the  stability  condition  for  the  wave  equation. 


At  < 


d-  +  \er,\  +  C^Tll  +  Tjl 


REFERENCES 


[1]  Anderson,  D.  A.,  Tannehill,  J.  C.,  and  Fletcher,  R.  H.,  Computational  Fluid 
Mechanics  and  Heat  Transfer,  Hemisphere  Publishing  Corporation,  New  York, 
1984. 

[2]  Anderson,  W.  Kyle,  Implicit  Multigrid  Algorithms  for  the  Three-Dimensional 
Flux  Split  Euler  Equations,  Ph.D.  Dissertation,  Mississippi  State  University, 
Starkville,  August  1986. 

[3]  Baysal,  0.,  Fouladi,  K.,  and  Lessard,  V.  R.,  “Multigrid  and  Upwind  Viscous  Flow 
Solver  on  Three-Dimensional  Overlapped  and  Embedded  Grids,”  AIAA  Journal, 
Vol.  29,  No.  6,  June  1991,  pp.  903-910. 

[4]  Belk,  Dave  M.,  Unsteady  Three-Dimensional  Euler  Equations  Solutions  on  Dy¬ 
namic  Blocked  Grids,  Ph.D.  Dissertation,  Mississippi  State  University,  Starkville, 
August  1986. 

[5]  Benek,  J.  A.,  Steger,  J.  L.,  Dougherty,  F.  C.,  and  Buning,  P.  G.,  “Chimera;  A 
Grid  Embedding  Technique,”  AEDC-TR-85-64,  December  1985. 

[6]  Berger,  M.  J.,  “On  Conservation  at  Grid  Interfaces,”  SIAM  Journal  on  Numerical 
Analysis,  Vol.  24,  No.  5,  October  1987,  pp.  967-984. 

[7]  Berger,  M.  J.  and  Jameson,  A.,  “Automatic  Adaptive  Grid  Refinement  for  the 
Euler  Equations,”  AIAA  Journal,  Vol.  23,  No.  4,  April  1985,  pp.  561-568. 

[8]  Brandt,  Achi,  “Multilevel  Adaptive  Computations  in  Fluid  Dynamics,”  AIAA 
Journal,  Vol.  18,  No.  10,  October  1980,  pp.  1165-1172. 

[9]  Brandt,  Achi,  “Multigrid  Techniques:  1984  Guide  with  Applications  to  Fluid 
Dynamics,”  Computational  Fluid  Dynamics  Lecture  Series  at  the  von-Karmen 
Institute  for  Fluid  Dynamics,  Rhode-Saint-Gense,  March  1984. 

[10]  Briggs,  William  L.,  A  Multigrid  Tutorial,  Society  for  Industrial  and  Applied 
Mathematics,  Philadelphia,  1987. 

[11]  Brown,  Jeffrey  J.,  “A  Multigrid  Mesh- Embedding  Technique  for  Three- 
Dimensional  Transonic  Potential  Flow  Analysis,”  NASA  Conference  Publication 
2202,  October  1981,  pp.  131-149. 

[12]  Dietz,  W.  E.  and  Suhs,  N.  E., “Pegasus  3.0  User’s  Manual,”  AEDC-TR-89-7, 
August  1989. 

[13]  Dougherty,  F.  Carroll,  Development  of  a  Chimera  Grid  Scheme  with  Applica¬ 
tions  to  Unsteady  Problems,  Ph.D.  Dissertation,  Stanford  University,  Palo  Alto, 
California,  June  1985. 

[14]  Dougherty,  F.  C.,  Benek,  J.  A.,  and  Steger,  J.  L.,  “On  Applications  of  Chimera 
Grid  Schemes  to  Store  Separation,”  NASA-TM-88193,  October  1985. 


70 


71 


[15]  Hackbusch,  Wolfgang,  Multi-Grid  Methods  and  Applications,  Springer- Verlag 
Berlin,  1985. 

[16]  Hanel,  D.,  Schroder,  W.,  and  Seider,  G.,  “Multigrid  Methods  for  the  Solution  of 
the  Compressible  Navier-Stokes  Equations,”  Proceedings  of  the  Fourth  GAMM- 
Seminar,  Kiel,  January  1988,  pp.  114-127. 

[17]  Henshaw,  W.  D.  and  Chesshirt,  G.,  “Multigrid  on  Composite  Meshes,”  SIAM 
Journal  on  Scientific  and  Statistical  Computation,  Vol.  8,  No.  6.  November  1987. 
pp.  914-923. 

[18]  Hessenius,  K.  A.  and  Rai,  M.  M.,  “Applications  of  a  Conservativ'^e  Zonal  Scheme 
to  Transient  and  Geometrically  Complex  Problems,”  AIAA-84-1532.  June  1984. 

[19j  John,  James  E.  A.,  Gas  Dynamics,  Second  Edition,  .A.llyn  and  Bacon,  Inc., 
Boston,  1984. 

[20]  Mavriplis,  D,  J.,  “Multigrid  Solution  of  the  Two  Dimensional  Euler  Equations 
on  Unstructured  Triangular  Meshes,”  AIAA  Journal,  Vol.  26.  No.  7.  Julv  1988, 
pp.  824-831. 

[21]  McCorm.ick,  Stephen  F.,  Multilevel  Adaptive  Methods  for  Partial  Differential 
Equations,  Society  for  Industrial  and  .Applied  Mathematics,  Philadelphia,  1989. 

[22]  Meakin,  R.  and  Suhs,  N.,  “Unsteady  .Aerodynamic  Simulation  of  Multiple  Bodies 
in  Relative  Motion,”  .AIAA-89-1996,  June  1989. 

[23]  Rai,  M.  M.,  “A  Conservative  Treatment  of  Zonal  Boundaries  for  Euler  Equation 
Calculations,”  AIAA-84-0164,  January  1984. 

[24]  Rai,  M.  M.,  “.An  Implicit,  Conservative,  Zonal-Boundary  Scheme  for  Euler  Equa¬ 
tion  Calculations,”  AIAA-85-0488,  January  1985. 

[25]  Roe,  P.  L.,  “.Approximate  Riemann  Solvers,  Parameter  Vectors,  and  Difference 
Schemes,”  Journal  of  Computational  Physics,  Vol.  43,  1981  pp.  357-372. 

[26]  Simpson,  L.  Bruce,  Unsteady  Three-Dimensional  Thin-Layer  Navier-Stokes  So¬ 
lutions  on  Dynamic  Blocked  Grids,  Ph.D.  Dissertation,  Mississippi  State  Univer¬ 
sity,  Starkville,  December  1988. 

[27]  Swanson,  R.  C.  and  Radespiel,  R.,  “Cell  Centered  and  Cell  Vertex  Multigrid 
Schemes  for  the  Navier-Stokes  Equations,”  .AIAA  Journal,  Vol.  29.  No.  5.  Mav 
1991,  pp.  697-703. 

[28]  Thompson,  J.  F.,  and  Gatlin,  B.,  “Program  E.AGLE  User's  Manual;  Surface 
Generation  Code,”  .AF.ATL-TR-88-117,  Vol.  II,  September  1988. 

[29]  Thompson,  J.  F.,  and  Gatlin,  B.,  “Program  EAGLE  User’s  Manual:  Grid  Gen¬ 
eration  Code,”  AFATL-TR-88-117,  Vol.  Ill,  September  1988. 

[30]  Thompson,  J.  F.,  Warsi,  Z.  .A.  U.,  and  Mastin,  C.  W.,  Numerical  Grid  Genera¬ 
tion:  Foundations  and  Applications,  North-Holland,  New  York,  1985. 

[31]  Whitfield,  D.  L.  and  Janus,  J.  M.,  “Three-Dimensional  Unsteady  Euler  Equations 
Solution  Using  Flux  Vector  Splitting,”  .AIAA-84-1552,  June  1984. 


72 


i32]  Whitfield,  D.  L.,  Janus,  J.  M.,  and  Simpson,  L.  B.,  “Implicit  Finite  Volume  High 
Resolution  Wave-Split  Scheme  for  Solving  the  Unsteady  Three-Dimensional  Euler 
and  Navier-Stokes  Equations  on  Stationary  or  Dynamic  Grids,”  MSSU-EIRS- 
ASE-88-2,  February  1988. 

33j  Yadlin,  Yoram  and  Caughey,  D.  A.,  “Block  Multigrid  Implicit  Solution  of  the 
Euler  Equations  of  Compressible  Fluid  Flow.”  AIAA  Journal.  Vol.  29.  No.  5. 
May  1991,  pp.  712-719. 


BIOGRAPHICAL  SKETCH 


Rudy  A.  Johnson  is  the  oldest  of  three  sons  born  to  Glen  and  Mary  Johnson. 
Born  in  Defiance,  Ohio,  in  1966  he  attended  Paulding  High  School  and  worked  part 
time  in  his  parents’  garage.  After  graduation  in  1984  he  attended  the  University  of 
Cincinnati  to  pursue  his  undergraduate  degree  in  aerospace  engineering.  In  1986  he 
began  working  as  a  co-op  student  in  the  computational  fluid  dynamics  (CFD)  section 
of  what  was  then  the  Air  Force  Armament  Laboratory  (AFATL/FXA)  at  Eglin  AFB. 
His  primary  responsibility  was  the  development  of  a  generalized  computer  graphics 
software  package  for  displaying  two  and  three  dimensional  CFD  results.  During  this 
time  he  also  co-authored  an  Air  Force  technical  report  (AFATL-TR-88-115). 

Upon  graduation  from  college  Ri  began  working  full  time  in  the  CFD  section  (cur¬ 
rently  designated  WL/MNAA)  as  a  junior  research  scientist /engineer-  and  program 
manager.  He  received  the  Junior  Scientist/Engineer  of  the  Year  Award”  in  May  1990 
and  has  received  performance  awards  each  year  since.  He  began  taking  graduate- 
level  classes  on  a  part  time  basis  during  the  fall  semester  of  1989  at  the  University  of 
Florida  Graduate  Research  Center  at  Eglin  AFB.  He  was  accepted  as  a  degree-seeking 
student  m  the  spring  of  1990.  In  October  of  1991  he  wed  Mary  Landers. 


7.3 


DISTRIBUTION 
(WL-TR -92-7039) 


Defense  Technical  Info.  Center 
Attn:  DTIC/DDAC  2 

Cameron  Station 
Alexandria  VA  22304-6145 

AFSAA/SAI  1 

Washington  DC  20330-5420 

AUL/LSE  1 

Maxwell  AFB  AL  36112-5564 


Eglin  AFB  offices: 

WUCA-N  1 

WL/MNOI  (Scientific  and  Tech.  Info.  Facility)  1 

HQ  USAFE/INATW  1 

APO  NY  0901 2-5001 

WUFIES/SURVIAC  1 

Wright-Fatterson  AFB  OH  45433-6553 


t 


Egiin  AFB  offices: 


WUMNSl  1 

WL/MNAG  1 

WUMNM  1 

ASC/XRC  1 

WUMNAA  4 

WUMNAV  1 

WUMNPX  1 


AFIA/INT  1 

Bolling  AFB  DC  20332-5000 

Commander 

U.  S.  Army  Missile  Command  ^ 

Redstone  Science  Info.  Center 
Attn:  AMSMI-RD-CS-R/Documents 
Redstone  Arsenal  AL  35898-5241 

Commander 

Naval  Weapons  Center  (Code  3431 )  ^ 

Attn:  Technical  Library 
China  Lake  CA  93555-6001 

NASA-Ames  Research  Center  ^ 

Attn:  Dr.  Tony  Strawa,  MS  229-3 
Moffett  Field  CA  94035-0000 


Wright-Patterson  AFB  OH  45433-6503 


ASC/ENSTA  1 

ASC/XRH  1 

Wright-Patterson  AFB  OH  45433-6553 

WUFIM  1 

WL7CA-F  1 

WL/FIB  1 

WLyFIGX  1 

WUFIGCC  1 

WL/TXA  1 

Wright-Patterson  AFB  OH  45433-6523 

EOARD/LDV  1 

Box  14 

FPO  NY  09510-0200 

NASA  Langley  Research  Center  1 

Technical  Library  -  MS  185 
Attn:  Document  Cataloging 
Hampton  VA  23665-5225 


74 


